'''
Replication Code for "How Students and the Public Define Terrorism,
and How Education Affects Those Definitions" 
'''



library(lme4)
library(tidyverse)
library(readxl)
library(writexl)
library(lavaan)
library(ggplot2)
library(semTools)
library(knitr)
library(semPlot)
library(stargazer)
library(lavaanPlot)
library(kableExtra)
library(officer)
library(flextable)
library(dplyr)


setwd("your own directory")
#Upload datasets
study1 <- read_excel( path = "Study1 dataverse.xlsx")
s2_exp<- read_excel(path = "Study2_experiment dataverse.xlsx")
s2_fup<-read_excel(path = "Study2_followup dataverse.xlsx")
study3 <- read_excel(path = "Study3 dataverse.xlsx")
study4 <- read_excel(path = "Study4 dataverse.xlsx")

## DESCRIPTIVE ----

sum(study1$time==0 & study1$treatment==1)
sum(study1$time==0 & study1$treatment==0)


sum(study3$time==0 & study3$treatment==1)
sum(study3$time==0 & study3$treatment==0)


sum(study4$time==0 & study4$treatment==1, na.rm=TRUE)
sum(study4$time==0 & study4$treatment==0, na.rm=TRUE)


median(as.numeric(substr(study4$Age, 1, 2)), na.rm=TRUE)
min(as.numeric(substr(study4$Age, 1, 2)), na.rm=TRUE)
max(as.numeric(substr(study4$Age, 1, 2)), na.rm=TRUE)

##### Percentage calculations  #####

variables_list = c("violence","state", "nsa", "targ_civ",  "targ_mil", "targ_pol","fear_terror", "obj_pol", 
                   "obj_any", "rel_any")


# Function to calculate percentage
calculate_percentage <- function(data, variable) {
  data %>%
    group_by(time, treatment) %>%
    summarise(percentage = mean(!!sym(variable), na.rm = TRUE) * 100,
              total_responses = sum(!is.na(!!sym(variable)))) %>%
    mutate(variable = variable)
}


percentage_s1 <- bind_rows(lapply(variables_list, calculate_percentage, data = study1)) %>%
  pivot_wider(names_from = variable, values_from = percentage) %>%
  mutate(Study = "Study 1")

# View the combined results
print(percentage_s1)

# study 2 is now study 3
percentage_s3 <- bind_rows(lapply(variables_list, calculate_percentage, data = study3)) %>%
  pivot_wider(names_from = variable, values_from = percentage) %>%
  mutate(Study = "Study 3")

percentage_s4 <- bind_rows(lapply(variables_list, calculate_percentage, data = study4)) %>%
  pivot_wider(names_from = variable, values_from = percentage) %>%
  mutate(Study = "Study 4")

calculate_percentage_s2 <- function(data, variable) {
  data %>%
    group_by(group_assignment) %>%
    summarise(percentage = mean(!!sym(variable), na.rm = TRUE) * 100,
              total_responses = sum(!is.na(!!sym(variable)))) %>%
    mutate(variable = variable)
}

## we have changed the names from study 4 to study 2
percentage_s2 <- bind_rows(lapply(variables_list, calculate_percentage_s2, data = s2_exp)) %>%
  pivot_wider(names_from = variable, values_from = percentage) %>%
  mutate(Study = "Study 2")

percentage_s2$treatment <- percentage_s2$group_assignment
perc_all <- rbind(percentage_s1, percentage_s2, percentage_s3, percentage_s4)


#### Academic and Government Definitions ####
getwd()
govt_media_acd <- read.csv("Government, Scholarship Definitions of Terrorism dataverse.csv", header = TRUE)
govt_media_acd[, 6:ncol(govt_media_acd)] <- lapply(govt_media_acd[, 6:ncol(govt_media_acd)], as.numeric)

calculate_percentage <- function(data, variable) {
  data %>%
    summarise(
      percentage = mean(!!sym(variable), na.rm = TRUE) * 100,
      total_responses = sum(!is.na(!!sym(variable)))
    ) %>%
    mutate(variable = variable)
}


study_profs <- govt_media_acd %>% 
  filter(Type == "Prof") 

study_profs_percentage <- bind_rows(lapply(variables_list, calculate_percentage, data = study_profs))
study_profs_percentage$Type <- "Study Profs"

# Filter for US govt and specific Source.Names
US_govt_def <- govt_media_acd %>%
  filter(Type == "US govt" & 
           Source.Name %in% c("U.S. State Department and CIA", 
                              "FBI (international terrorism)", 
                              "US Department of Homeland Security"))

usandfor_govt_def <- govt_media_acd %>%
  filter(Type == "US govt and Foreign govt")

US_govt_def <- bind_rows(US_govt_def, usandfor_govt_def)
US_govt_percentage <- bind_rows(lapply(variables_list, calculate_percentage, data = US_govt_def))
US_govt_percentage$Type <- "US Government"


# Filter for Foreign govt, Intl org, and US govt and Foreign govt without Source.Name restriction
nonUS_govt_def <- govt_media_acd %>%
  filter(Type %in% c("Foreign govt", "Intl org"))

nonUS_percentage <- bind_rows(lapply(variables_list, calculate_percentage, data = nonUS_govt_def))
nonUS_percentage$Type <- "non-US Government"


# Combine the two datasets
govt_def <- bind_rows(US_govt_def, nonUS_govt_def)
govt_percentage <- bind_rows(lapply(variables_list, calculate_percentage, data = govt_def))
govt_percentage$Type <- "Government"

academic_def <- govt_media_acd %>% 
  filter(Type == "Academic" |Type == "Academic (military)" ) %>% 
  drop_na(all_of(variables_list))

academic_percentage <- bind_rows(lapply(variables_list, calculate_percentage, data = academic_def))
academic_percentage$Type <- "Academic"

# Filter for US govt and specific Source.Names
US_govt_pre911_def <- govt_media_acd %>%
  filter(Type == "US govt" & Year < 2001)
US_govt_pre911_percentage <- bind_rows(lapply(variables_list, calculate_percentage, data = US_govt_pre911_def))
US_govt_pre911_percentage$Type <- "US government pre 9/11"

Perc_def <- rbind(govt_percentage,US_govt_percentage, nonUS_percentage, academic_percentage, study_profs_percentage, US_govt_pre911_percentage)

# Save the new data frame to a CSV file if needed
#write_xlsx(Perc_def, "media_academic_govt_perc.xlsx")

## Study 1 Difference in Difference Analysis  2. Regression -----

dependent_variables<- c("violence", "state", "nsa",  "targ_civ", "targ_pol", "targ_mil",  "fear_terror", "obj_pol", "obj_any", "rel_any", "rel_isl" )

models <- list()

# Run regression for each dependent variable and store models in the list

for (dep_var in dependent_variables) {
  model <- lm(paste(dep_var, "~ time * treatment + knowledge + interest "), data = study1)
  models[[dep_var]] <- model
}


table<- stargazer(models, type = "html", digits = 2, title = "Estimates of difference-in-difference analysis (Study 1)",
                  column.labels = c("Violence", "State", "NSA",  "Civilian target", "Political target", "Military target", "Fear or Terror", "Political objective",
                                    "Any objectives", "Religion", "Islam"),
                  dep.var.caption  = "Dependent Variables",
                  covariate.labels=c("time", "treatment", "knowledge", "interest", "time*treatment", "Constant"), 
                  out = "s1_regression_table.html")
#print(table)


## Study 1 Difference in Difference Analysis  1. Plot -----


# Create a data frame to store the results
results <- data.frame(
  Dependent_Variable = character(),
  Estimate = numeric(),
  Std_Error = numeric(),
  P_Value = numeric()
)

# Loop over each dependent variable, run DiD regression, and store the estimate, standard error, and p-value
for (variable in dependent_variables) {
  # Run DiD regression
  did_model <- lm(paste(variable, "~ time * treatment + knowledge + interest"), data = study1)
  
  # Extract coefficient for the treatment effect (interaction term)
  coefficient <- coef(did_model)["time:treatment"]
  std_error <- sqrt(vcov(did_model)["time:treatment", "time:treatment"])  # Standard error
  p_value <- summary(did_model)$coefficients["time:treatment", "Pr(>|t|)"]  # P-value
  
  # Store results in the data frame
  results <- rbind(results, data.frame(Dependent_Variable = variable,
                                       Estimate = coefficient,
                                       Std_Error = std_error,
                                       P_Value = p_value))
}

results$Dependent_Variable <- factor(results$Dependent_Variable, levels = reversed_order)
#filtered_results <- results %>%
#  filter(!Dependent_Variable %in% excluded_variables)

# Create the plot
did_plot_s1 <- ggplot(results[!is.na(results$Estimate), ], aes(x = Estimate, y =Dependent_Variable)) +
  geom_point(aes(color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), size = 3) +
  geom_line(aes(group = Dependent_Variable, color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), size = 1) +
  geom_errorbarh(aes(xmax = Estimate + Std_Error, xmin = Estimate - Std_Error, color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), height = 0) +  # Change error bar color based on significance
  labs(title = "Study 1",
       x = "Estimate",
       y = "Dependent Variable") +
  scale_color_manual(values = c("Significant" = "black", "Non-Significant" = "gray")) +  # Set colors for significant and non-significant
  theme_minimal()+
  theme(legend.position = "none")+
  theme(plot.background = element_rect(fill = "white"))

# Display the plot
print(did_plot_s1)
#ggsave("s1_did_plot.png", did_plot, width = 8, height = 6, dpi = 300)

## Study 4 Difference in Difference Analysis  1. Plot ####

for (variable in dependent_variables) {
  if (variable %in% colnames(study4)) {  # Check if the variable exists in the dataset
    # Try converting the variable to numeric
    study4[[variable]] <- as.numeric(study4[[variable]])
    
    # Check for NAs after conversion
    if (any(is.na(study4[[variable]]))) {
      cat("Warning: Variable", variable, "introduced NAs after coercion.\n")
      cat("Number of NAs introduced:", sum(is.na(study4[[variable]])), "\n")
    }
  } else {
    cat("Variable", variable, "not found in the dataset.\n")
  }
}

results <- data.frame(
  Dependent_Variable = character(),
  Estimate = numeric(),
  Std_Error = numeric(),
  P_Value = numeric()
)

# Loop over each dependent variable, run DiD regression, and store the estimate, standard error, and p-value

for (variable in dependent_variables) {
  # Check if the variable exists in the dataset
  if (variable %in% colnames(study4)) {
    
    # Subset data to remove rows with NA in the dependent variable and the other model variables
    temp_data <- study4[!is.na(study4[[variable]]) &
                          !is.na(study4$time) &
                          !is.na(study4$treatment) &
                          !is.na(study4$knowledge) &
                          !is.na(study4$interest), ]
    
    # Run DiD regression if there is enough data left after removing NAs
    if (nrow(temp_data) > 0) {
      did_model <- lm(paste(variable, "~ time * treatment + knowledge + interest"), data = temp_data)
      
      # Extract coefficient for the treatment effect (interaction term)
      coefficient <- coef(did_model)["time:treatment"]
      std_error <- sqrt(vcov(did_model)["time:treatment", "time:treatment"])  # Standard error
      p_value <- summary(did_model)$coefficients["time:treatment", "Pr(>|t|)"]  # P-value
      
      # Store results in the data frame
      results <- rbind(results, data.frame(Dependent_Variable = variable,
                                           Estimate = coefficient,
                                           Std_Error = std_error,
                                           P_Value = p_value))
    } else {
      cat("Not enough data for variable", variable, "after removing NAs.\n")
    }
  } else {
    cat("Variable", variable, "not found in the dataset.\n")
  }
}



desired_order <- c("violence","state", "nsa","targ_civ", "targ_pol", "targ_mil",  "fear_terror", "fear_end", "fear_means",
                   "obj_pol", "obj_ideo", "obj_any", "rel_any", "rel_isl")
reversed_order <- rev(desired_order)
results$Dependent_Variable <- factor(results$Dependent_Variable, levels = reversed_order)

# Exclude the unwanted variables
#excluded_variables <- c("fear_end", "fear_means","obj_ideo")

# Filter the results to exclude the unwanted variables
#filtered_results <- results %>%
#  filter(!Dependent_Variable %in% excluded_variables)

# Create the plot with the filtered data
did_plot_s4 <- ggplot(results[!is.na(results$Estimate), ], aes(x = Estimate, y = Dependent_Variable)) +
  geom_point(aes(color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), size = 3) +
  geom_line(aes(group = Dependent_Variable, color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), size = 1) +
  geom_errorbarh(aes(xmax = Estimate + Std_Error, xmin = Estimate - Std_Error, color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), height = 0) +
  labs(title = "Study 4",
       x = "Estimate",
       y = "Dependent Variable") +
  scale_color_manual(values = c("Significant" = "black", "Non-Significant" = "gray")) +
  theme_minimal() +
  theme(legend.position = "none") +
  theme(plot.background = element_rect(fill = "white"))
# Display the plot
print(did_plot_s4)

##DiD Plots  ####                      
#Dependent variables 
dependent_variables <- c("violence","state", "nsa","targ_civ", "targ_pol", "targ_mil",  "fear_terror",
                         "obj_pol", "obj_any", "rel_any", "rel_isl")

# Define the desired order of dependent variables and reverse it
desired_order <- c('rel_isl', 'rel_any', 'obj_any', 'obj_pol', 
                   'fear_terror', 'targ_mil', 'targ_pol', 'targ_civ', 'nsa', 
                   'state', 'violence')

# Create a named vector that maps variable names to labels
variable_labels <- c(
  "violence" = "Violence",
  "state" = "State",
  "nsa" = "Non-State Actor",
  "targ_civ" = "Civilian Target",
  "targ_pol" = "Political Target",
  "targ_mil" = "Military Target",
  "fear_terror" = "Fear or Terror",
  "obj_pol" = "Political Objective",
  "obj_any" = "Any Objective",
  "rel_any" = "Religion Any",
  "rel_isl" = "Religion Islam"
)

# Create a data frame to store the results
results <- data.frame(
  Dependent_Variable = character(),
  Estimate = numeric(),
  Std_Error = numeric(),
  P_Value = numeric()
)

# Function to store DiD results for each study
store_did_results <- function(data, study_name, model_formula, desired_order) {
  results <- data.frame(
    Study = character(),
    Dependent_Variable = character(),
    Estimate = numeric(),
    Std_Error = numeric(),
    P_Value = numeric()
  )
  
  for (variable in dependent_variables) {
    # Run DiD regression
    did_model <- lm(as.formula(paste(variable, model_formula)), data = data)
    coefficient <- coef(did_model)["time:treatment"]
    std_error <- sqrt(vcov(did_model)["time:treatment", "time:treatment"])
    p_value <- summary(did_model)$coefficients["time:treatment", "Pr(>|t|)"]
    
    results <- rbind(results, data.frame(
      Study = study_name,
      Dependent_Variable = variable,
      Estimate = coefficient,
      Std_Error = std_error,
      P_Value = p_value
    ))
  }
  
  # Ensure the Dependent_Variable column is a factor with the specified order
  results$Dependent_Variable <- factor(results$Dependent_Variable, levels = desired_order)
  return(results)
}


# Store results for each study
results_s1 <- store_did_results(study1, "Study 1", "~ time * treatment + knowledge + interest", desired_order)
results_s2 <- store_did_results(study3, "Study 3", "~ time * treatment +  knowledge + interest", desired_order)
results_s3 <- store_did_results(study4, "Study 4", "~ time * treatment  +  knowledge + interest", desired_order)
combined_results <- bind_rows(results_s1, results_s2, results_s3)
combined_results$Dependent_Variable <- recode(combined_results$Dependent_Variable, !!!variable_labels)



# Create the plot with facet_wrap
# Create the plot with larger fonts and custom labels
combined_plot <- ggplot(combined_results[!is.na(combined_results$Estimate), ], aes(x = Estimate, y = Dependent_Variable)) +
  geom_point(aes(color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), size = 3) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  geom_errorbarh(aes(xmax = Estimate + Std_Error, xmin = Estimate - Std_Error, color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), height = 0) +
  labs(title = "Estimates of Difference-in-Differences Analysis",
       x = "Estimate",
       y = "Dependent Variable") +
  scale_color_manual(values = c("Significant" = "black", "Non-Significant" = "gray")) +
  theme_minimal() +
  theme(
    legend.position = "none",
    plot.background = element_rect(fill = "white"),
    plot.title = element_text(size = 14, face = "bold"),           # Increase title font size
    axis.title = element_text(size = 14),                         # Increase axis titles font size
    axis.text = element_text(size = 14),                          # Increase axis text font size
    strip.text = element_text(size = 15, face = "bold"),                         # Increase facet labels font size
    legend.text = element_text(size = 14)                         # Increase legend text font size (if visible)
  ) +
  facet_wrap(~ Study, scales = "free_x")

# Display the plot
print(combined_plot)

#ggsave("s1_3_4_did_plots.png", combined_plot, width = 8, height = 8, dpi = 600)


## Study 4 LINEAR REGRESSION  1. Plot  --------- 

# Define the function to store DiD results
store_lm_results <- function(data, study_name) {
  results <- data.frame(
    Dependent_Variable = character(),
    Estimate = numeric(),
    Std_Error = numeric(),
    P_Value = numeric(),
    Study = character()
  )
  
  dependent_variables <-  c( "violence", "state", "nsa", "fear_terror","targ_civ", "targ_pol", "obj_pol", "obj_any", "rel_any")
  
  for (variable in dependent_variables) {
    # Run DiD regression
    lm_model <- lm(as.formula(paste(variable, "~ group_assignment")), data = data)
    
    # Extract coefficient for the treatment effect (interaction term)
    coefficient <- coef(lm_model)["group_assignment"]
    std_error <- sqrt(vcov(lm_model)["group_assignment", "group_assignment"])  # Standard error
    p_value <- summary(lm_model)$coefficients["group_assignment", "Pr(>|t|)"]  # P-value
    
    # Store results in the data frame
    results <- rbind(results, data.frame(
      Dependent_Variable = variable,
      Estimate = coefficient,
      Std_Error = std_error,
      P_Value = p_value,
      Study = study_name
    ))
  }
  
  return(results)
}

# Store results for experiment and follow-up data
results_exp <- store_lm_results(s2_exp, "Experiment")
results_followup <- store_lm_results(s2_fup, "Follow-Up")

# Combine the results
s4_results <- bind_rows(results_exp, results_followup)


desired_order <- c( "violence", "state", "nsa", "targ_civ", "targ_pol", "fear_terror", "obj_pol", "obj_any", "rel_any")
reversed_order <- rev(desired_order)
s4_results$Dependent_Variable <- factor(s4_results$Dependent_Variable, levels = reversed_order)
s4_results$Dependent_Variable <- recode(s4_results$Dependent_Variable, !!!variable_labels)

# Create the s4 plot
s4_plot <- ggplot(s4_results[!is.na(s4_results$Estimate), ], aes(x = Estimate, y = Dependent_Variable)) +
  geom_point(aes(color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), size = 3) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  geom_errorbarh(aes(xmax = Estimate + Std_Error, xmin = Estimate - Std_Error, color = ifelse(P_Value < 0.05, "Significant", "Non-Significant")), height = 0) +
  labs(title = "Estimates of Linear Regression Analysis (Study 2)",
       x = "Estimate",
       y = "Dependent Variable") +
  scale_color_manual(values = c("Significant" = "black", "Non-Significant" = "gray")) +
  theme_minimal() +
  theme(
    legend.position = "none",
    plot.background = element_rect(fill = "white"),
    plot.title = element_text(size = 14, face = "bold"),           # Increase title font size
    axis.title = element_text(size = 14),                         # Increase axis titles font size
    axis.text = element_text(size = 14),                          # Increase axis text font size
    strip.text = element_text(size = 15, face = "bold"),                         # Increase facet labels font size
    legend.text = element_text(size = 14)                         # Increase legend text font size (if visible)
  ) +
  theme(legend.position = "none") +
  theme(plot.background = element_rect(fill = "white")) +
  facet_wrap(~ Study, scales = "free_x")

# Display the plot
print(s4_plot)


## Latent variable analysis ------

model <- '
  # Measurement model for the Rational variable
  Rational =~   fear_means + obj_pol +  targ_pol + obj_ideo + targ_mil + obj_org

  # Measurement model for the Irrational variable
  Irrational =~ fear_end

  # Structural model to test the impact of the experiment
  Rational ~ c * time + treatment + did 
  Irrational ~ d *  time + treatment + did 

'

# Fit the model
fit_study1 <- cfa(model, data = study1)

# Create the path diagram
#semPaths(fit_study1, nCharNodes = 10)

fit_study3 <- cfa(model, data = study3)

# Create the path diagram
#semPaths(fit_study3, nCharNodes = 10)

fit_study4 <- cfa(model, data = study4)

# Create the path diagram
#semPaths(fit_study4, nCharNodes = 10)

summary(fit_study1)
summary(fit_study3)
summary(fit_study4)

model_s2 <- '
  # Measurement model for the Rational variable
  Rational =~ fear_means + obj_pol +  targ_pol + obj_ideo + targ_mil +obj_org

  # Measurement model for the Irrational variable
  Irrational =~ fear_end

  # Structural model to test the impact of the experiment
  Rational ~ c * group_assignment 
  Irrational ~ d * group_assignment
'

colnames(s2_exp)

# Fit the model
fit_study2 <- cfa(model_s2, data = s2_exp)
fit_study2_fup <- cfa(model_s2, data = s2_fup)

summary(fit_study2)

labels <- NULL
# Create the path diagram
labels$fear_means<- "Fear as Means"
labels$obj_pol <- "Political Objective"
labels$obj_ideo <- "Ideological Objective"
labels$targ_mil <- "Military Target"
labels$obj_ideo <- "Ideological Objective"
labels$obj_org <- "Organizational Objective"
labels$targ_pol <- "Political Target"
labels$did <- "Treatment * Time"
labels$fear_end <- "Fear as an End"
labels$group_assignment <-NULL

lavaanPlot(fit_study1,  labels = labels)

labels$did <-NULL
labels$group_assignment <-"Treatment"
lavaanPlot(fit_study2,  labels = labels)

# study 1 Lavaan findings:
labels$did <- "Time*Treatment"


s1_lp_circo<-lavaanPlot(model = fit_study1, coefs = TRUE, labels = labels,  stars = "regress",
                        node_options = list(shape = "box", fontname = "TimesNewRoman"), 
                        edge_options = list(color = "black")) 


print(s1_lp_circo)
# Extract standardized factor loadings
factor_loadings <- parameterEstimates(fit_study1, standardized = TRUE)
factor_loadings <- subset(factor_loadings, op == "=~")

# factor loadings (correlations of indicators with latent variables)
fit_table_html <- kable(
  factor_loadings[, c("lhs", "rhs", "est", "se", "pvalue")], 
  format = "html", 
  caption = "Factor Loadings and Model Parameters - Study 1", 
  col.names = c("Latent Factor", "Observed Variable","Parameter Value", "Std. Error", "p-Value")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

# Save the HTML table to a file
save_kable(fit_table_html, file = "Latent_fit_study1.html")

# Create a flextable object from the data frame
fit_flextable <- flextable(factor_loadings[, c("lhs", "rhs", "est", "se", "pvalue")])
fit_flextable <- set_header_labels(fit_flextable, 
                                   lhs = "Latent Factor", 
                                   rhs = "Observed Variable", 
                                   est = "Parameter Value", 
                                   se = "Std. Error", 
                                   pvalue = "p-Value")

fit_flextable <- autofit(fit_flextable)

# Create a new Word document
doc <- read_docx()

# Add title and the table to the document
doc <- doc %>%
  body_add_par("Factor Loadings and Model Parameters - Study 1", style = "heading 1") %>%
  body_add_flextable(fit_flextable)

# Save the Word document
print(doc, target = "fit_study1.docx")


# study 3 (university sample) Lavaan findings:

s3_lp_circo<- lavaanPlot(model = fit_study3, coefs = TRUE, labels = labels,  stars = "regress",
                         node_options = list(shape = "box", fontname = "TimesNewRoman"), 
                         edge_options = list(color = "darkgrey")) 

print(s3_lp_circo)

factor_loadings <- parameterEstimates(fit_study3, standardized = TRUE)
factor_loadings <- subset(factor_loadings, op == "=~")

# factor loadings (correlations of indicators with latent variables)
fit_table_html <- kable(
  factor_loadings[, c("lhs", "rhs", "est", "se", "pvalue")], 
  format = "html", 
  caption = "Factor Loadings and Model Parameters - Study 2", 
  col.names = c("Latent Factor", "Observed Variable","Parameter Value", "Std. Error", "p-Value")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

# Save the HTML table to a file
save_kable(fit_table_html, file = "Latent_fit_study3.html")

# Create a flextable object from the data frame
fit_flextable <- flextable(factor_loadings[, c("lhs", "rhs", "est", "se", "pvalue")])
fit_flextable <- set_header_labels(fit_flextable, 
                                   lhs = "Latent Factor", 
                                   rhs = "Observed Variable", 
                                   est = "Parameter Value", 
                                   se = "Std. Error", 
                                   pvalue = "p-Value")

fit_flextable <- autofit(fit_flextable)

# Create a new Word document
doc <- read_docx()

# Add title and the table to the document
doc <- doc %>%
  body_add_par("Factor Loadings and Model Parameters - Study 3", style = "heading 1") %>%
  body_add_flextable(fit_flextable)

# Save the Word document
print(doc, target = "fit_study4(uni).docx")

# study 4 findings
s4_lp_circo<- lavaanPlot(model = fit_study4, coefs = TRUE, labels = labels,  stars = "regress",
                         node_options = list(shape = "box", fontname = "TimesNewRoman"), 
                         edge_options = list(color = "darkgrey")) 
print(s4_lp_circo)

factor_loadings <- parameterEstimates(fit_study4, standardized = TRUE)
factor_loadings <- subset(factor_loadings, op == "=~")

# factor loadings (correlations of indicators with latent variables)
fit_table_html <- kable(
  factor_loadings[, c("lhs", "rhs", "est", "se", "pvalue")], 
  format = "html", 
  caption = "Factor Loadings and Model Parameters - Study 3", 
  col.names = c("Latent Factor", "Observed Variable","Parameter Value", "Std. Error", "p-Value")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

# Save the HTML table to a file
save_kable(fit_table_html, file = "Latent_fit_study4.html")


# Create a flextable object from the data frame
fit_flextable <- flextable(factor_loadings[, c("lhs", "rhs", "est", "se", "pvalue")])
fit_flextable <- set_header_labels(fit_flextable, 
                                   lhs = "Latent Factor", 
                                   rhs = "Observed Variable", 
                                   est = "Parameter Value", 
                                   se = "Std. Error", 
                                   pvalue = "p-Value")

fit_flextable <- autofit(fit_flextable)

# Create a new Word document
doc <- read_docx()

# Add title and the table to the document
doc <- doc %>%
  body_add_par("Factor Loadings and Model Parameters - Study 3", style = "heading 1") %>%
  body_add_flextable(fit_flextable)

# Save the Word document
print(doc, target = "fit_study2(mooc).docx")

# study 2 findings 
labels$did <- NULL
labels$group_assignment <- "Treatment"

s2_lavaanPlot<-lavaanPlot(model = fit_study2, coefs = TRUE, labels = labels,  stars = "regress",
                          node_options = list(shape = "box", fontname = "TimesNewRoman"), 
                          edge_options = list(color = "black")) 

print(s2_lavaanPlot)

s2_lavaanPlot<-lavaanPlot(model = fit_study2_fup, coefs = TRUE, labels = labels,  stars = "regress",
                          node_options = list(shape = "box", fontname = "TimesNewRoman"), 
                          edge_options = list(color = "grey")) 

print(s2_lavaanPlot)


factor_loadings <- parameterEstimates(fit_study2, standardized = TRUE)
factor_loadings <- subset(factor_loadings, op == "=~")

# factor loadings (correlations of indicators with latent variables)
fit_table_html <- kable(
  factor_loadings[, c("lhs", "rhs", "est", "se", "pvalue")], 
  format = "html", 
  caption = "Factor Loadings and Model Parameters - Study 2", 
  col.names = c("Latent Factor", "Observed Variable","Parameter Value", "Std. Error", "p-Value")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

# Save the HTML table to a file
save_kable(fit_table_html, file = "Latent_fit_study3_mturk.html")


# Create a flextable object from the data frame
fit_flextable <- flextable(factor_loadings[, c("lhs", "rhs", "est", "se", "pvalue")])
fit_flextable <- set_header_labels(fit_flextable, 
                                   lhs = "Latent Factor", 
                                   rhs = "Observed Variable", 
                                   est = "Parameter Value", 
                                   se = "Std. Error", 
                                   pvalue = "p-Value")

fit_flextable <- autofit(fit_flextable)

# Create a new Word document
doc <- read_docx()

# Add title and the table to the document
doc <- doc %>%
  body_add_par("Factor Loadings and Model Parameters - Study 2", style = "heading 1") %>%
  body_add_flextable(fit_flextable)

# Save the Word document
print(doc, target = "fit_table.docx")




#### Model Fit ####
model3<- '
  # Measurement model for the Rational variable
  Rational =~ fear_means + obj_pol +  targ_pol + obj_ideo + targ_mil +obj_org

  # Measurement model for the Irrational variable
  Irrational =~ fear_end

  # Structural model to test the impact of the experiment
  Rational ~ c * group_assignment 
  Irrational ~ d * group_assignment
'

model2 <- '
  # Measurement model for the Rational variable
  Rational =~ fear_means + obj_pol + obj_ideo+  targ_pol

  # Measurement model for the Irrational variable
  Irrational =~ fear_end

  # Structural model to test the impact of the experiment
  Rational ~ c * group_assignment 
  Irrational ~ d * group_assignment
'

model1 <- '
  # Measurement model for the Rational variable
  Rational =~ fear_means + obj_pol +  targ_pol

  # Measurement model for the Irrational variable
  Irrational =~ fear_end

  # Structural model to test the impact of the experiment
  Rational ~ c * group_assignment 
  Irrational ~ d * group_assignment
'

fıt_s2_1 <- cfa(model1, data = s2_exp)
fıt_s2_2 <- cfa(model2, data = s2_exp)
fıt_s2_3 <- cfa(model3, data = s2_exp)

# Load necessary libraries
library(semTools)

latent_variable_estimates <- inspect(fit_study3, "start")
standard_errors <- inspect(fit_study3, "se")

# Perform likelihood ratio test
lrt_result <- lavTestLRT(fıt_s2_1, fıt_s2_2, fıt_s2_3)
# Create a table
rownames(lrt_result) <- c("Model 1","Model 2", "Model 3")
lrt_table <- kable(lrt_result, format = "html", caption = "Likelihood Ratio Test Results") %>%
  kable_styling(full_width = F)

# Save to HTML
writeLines(lrt_table, "lrt_results.html")

# Define function to extract fit indices
extract_fit_indices <- function(fit) {
  fitMeasures(fit, c("chisq", "df", "pvalue", "cfi", "rmsea", "aic", "bic"))
}

# Extract fit indices for each model
fit_indices <- lapply(list(fıt_s2_1, fıt_s2_2, fıt_s2_3), extract_fit_indices)

# Create a data frame to store model fit information
fit_table <- data.frame(
  Model = c("Model 1", "Model 2", "Model 3"),
  ChiSquare = sapply(fit_indices, function(x) x["chisq"]),
  df = sapply(fit_indices, function(x) x["df"]),
  pvalue = sapply(fit_indices, function(x) x["pvalue"]),
  CFI = sapply(fit_indices, function(x) x["cfi"]),
  RMSEA = sapply(fit_indices, function(x) x["rmsea"]),
  AIC = sapply(fit_indices, function(x) x["aic"]),
  BIC = sapply(fit_indices, function(x) x["bic"])
)

# View the fit table
print(fit_table)

# Create a flextable object from the data frame
fit_flextable <- flextable(fit_table)
fit_flextable <- autofit(fit_flextable)
# Create a new Word document
doc <- read_docx()

# Add title and the table to the document
doc <- doc %>%
  body_add_par("Model Fit Summary", style = "heading 1") %>%
  body_add_flextable(fit_flextable)

# Save the Word document
print(doc, target = "fitmodel_table.docx")


# Print the model fit table
kable(fit_table, caption = "Model Fit Summary", out = "Model_fit.html")


# Create the table with kable and customize it with kableExtra
fit_table_html <- kable(
  fit_table, 
  format = "html", 
  caption = "Model Fit Summary", 
  col.names = c("Model", "Chi-Square", "df", "p-value", "CFI", "RMSEA", "AIC", "BIC")
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

# Save the HTML table to a file
save_kable(fit_table_html, file = "Model_fit.html")


print(modindices(fit_study1))


